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In this paper, we present an efficient and stable method to determine the one-particle Green's 
function in the hybridization-expansion continuous-time (CT-HYB) quantum Monte Carlo method, 
within the framework of the dynamical mean-field theory (DMFT). The high-frequency tail of the 
impurity self-energy is replaced with a noise-free function determined by a dual-expansion around 
the atomic limit. This method does not depend on the explicit form of the interaction term. More 
advantageous, it does not introduce any additional numerical cost to the CT-HYB simulation. We 
discuss the symmetries of the two-particle vertex, which can be used to optimize the simulation of the 
four-point correlation functions in the CT-HYB. Here, we adopt it to accelerate the dual-expansion 
calculation, which turns out to be especially suitable for the study of material systems with com- 
plicated band structures. As an application, a two-orbital Anderson impurity model with a general 
on-site interaction form is studied. The phase diagram is extracted as a function of the Coulomb 
interactions for two different Hund's coupling strengths. In the presence of the hybridization be- 
tween different orbitals, for smaller interaction strengths, this model shows a transition from metal 
to band-insulator. Increasing the interaction strengths, this transition is replaced by a crossover 
from Mott-insulator to band-insulator behavior. 

PACS numbers: 71.10.Fd, 71.27. +a, 71.30.+h 



I. INTRODUCTION 

The study of electronic structure of transition metal 
and heavy-fermion materials is one of the most active 
fields in condensed-matter physics. The highly correlated 
d- and /-electrons cannot be fully described by effective 
single-particle methods, such as the local-density approx- 
imation (LDA) to the density- functional theory (DFT). 
Here, the dynamical mean-field theory (DMFT) can be 
a powerful tool, especially when the momentum depen- 
dence of the self-energy is essentially negligible, regard- 
less of the electron-electron interaction strength. 1-3 

The central problem of the DMFT is to solve an effec- 
tive impurity model. In real materials, such a model usu- 
ally contains both inter- and intra-orbital interactions, 
as well as the hybridization among different orbitals. 
They account for the competitions between the mag- 
netic, charge, and orbital fluctuations. Thus, an efficient 
impurity solver, which can handle all the interactions 
and hybridizations, is of obvious importance. Among 
the available impurity solvers, 4-8 the numerically exact 
quantum Monte Carlo (QMC) methods were widely used. 
The recent development of the continuous-time quantum 
Monte Carlo (CT-QMC) methods 9 " 12 further supports 
the DMFT for the study of realistic materials in the sense 
that lower temperature regions can be reached and more 
orbitals can be investigated. 

For realistic material calculations based on the CT- 
QMC solvers, correctly resolving the high-frequency be- 
havior of the impurity self-energy T, imp (iu; n ) is of crucial 
importance. On the one hand, due to the iterative na- 
ture of the DMFT equations, T,i mp (iuj m ) determines the 
Weiss function at each iteration and, in the end, the con- 
verged solution of the DMFT procedure in some cases. 



On the other hand, T,i mp (iuj n ) strongly influences the de- 
termination of the total particle number and the analyt- 
ical continuation for a full spectral function calculation, 
which has a direct connection to experiments. 

In this paper, we show how to determine the impu- 
rity self-energy for a rather general multiorbital model in 
an efficient and stable manner within the hybridization- 
expansion continuous-time (CT-HYB) method. The di- 
rect simulation in the Matsubara-frequency space and 
careful treatment of the self-energy high-frequency tail 
make this method especially suitable for studying the 
material systems with complex band structures. 

This paper is organized as follow: Sec. II explains how 
the "dual transformation" can be employed to simulate 
effectively the one-particle Green's function in the CT- 
HYB. Additionally, it is shown how the simulation of the 
two-particle Green's function \ can be straightforwardly 
carried out as Wick's theorem still holds. The symmetry 
of x is discussed in detail in this section. In Sec. Ill, we 
make use of the CT-HYB to study a two-orbital Hubbard 
model with a general interaction term. For readers who 
are especially interested in our CT-HYB implementation 
and the self-energy correction scheme, Sec. II is the pri- 
mary option. If the phase diagram of the two-orbital 
model is of primary interest, one may skip Sec. II and 
go to Sec. Ill, which is self-contained. Conclusions and 
outlook can be found in Sec. IV. 



II. METHOD 

To explain our implementation of the CT-HYB in a 
concrete framework, we take a two-orbital model as an 
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example, that is, 

Hloc = H A + Hy + H int — 22 (1) 

0=1,2 a 

where H& = A ^2 a {ni r7 — n2 a ) represents the crystal field 
splitting and Hy = V X^o-( c i(T C 2cr + h.c.) is the hybridiza- 
tion between two orbitals. For the interaction part, a 
general on-site form is considered, 



-7(4^4^024.01-)- + c^c^ci t c u + h.c.) 



t J 



(2) 



which contains the intraorbital and interorbital Coulomb 
interactions, as well as the spin-flip and pair-hopping pro- 
cesses. 

As an impurity solver for the DMFT, CT-HYB cm- 
ploys the same idea as all the other CT-QMC impurity 
solvers; that is, it expands the impurity effective action 
around a certain limit and evaluates the expansion terms 
via stochastic sampling. Here, we only present the ex- 
pressions relevant to this work. For a more detailed re- 
view of the CT-QMC methods, we suggest Ref. 13. 

In the CT-HYB, the expansion of the "impurity + 
bath" action Stot — Si oc + Sbath + Shyb around the atomic 
limit is carried out by integrating out the bath degrees of 
freedom. Si oc , Sbath are the actions for the local and the 
bath Hamiltonian, respectively. Sh y b is the hybridization 
between them, 3 which is expanded order by order. The 
contraction of the bath operator b a , b\ follows Wick's the- 
orem, as the bath is noninteracting. This results in a de- 
terminant Det Ck with the hybridization function A(r, r') 
as matrix elements. The full hybridization matrix usually 
can be decoupled into block diagonal form with respect 
to certain conserved quantum numbers, for example, the 
total particle number n, the spin a z and cluster momenta 
K. The final expression of the partition function can then 
be written as 



ZbZioc TT J2 / [1 dT * dT * MCk a )Det c 

a k a J ° i=l 



(3) 



Here, k a is the expansion order (also the dimension of 
the determinant matrix) for the "a" flavor, where flavor 
represents spin, orbital, or cluster momenta. Tr(C^) = 
( T r Ua c a( T i)4( r i) ' ' ' c aK)4( T fc)) is the cluster trace of 
a group of "kinks", 12 that is, cluster operators, in the in- 
terval [0, j3). From now on, we always work with the diag- 
onal form of the hybridization function. The evaluation 
of Tr(Ck a ) can be carried out in two ways. One can either 
express the c a , c\ operators as matrices in the eigenba- 
sis of Hioc or employ the Krylov implementation 14 . The 
former one benefits from the diagonal form of the time 
evolution operators e~ HlocT . The Krylov implementa- 
tion, on the other hand, works in the particle-number 
basis, for which e~ Hl ° cT becomes a sparse matrix. It uses 
the efficient Krylov-space method, which makes it pos- 
sible to simulate up to typically seven orbital problems 



at acceptable numerical costs. In this work, the first im- 
plementation is used, in which we diagonalize Hi oc with 
respect to the conserved quantum numbers. 15 The trace 
of the fermion operators is evaluated by first searching 
for nonzero overlap between different eigenstates with re- 
spect to the group of the cluster operators. The nonzero 
trace is, then, calculated along the trajectory found. 



A. One-particle Green's function 

The impurity Green's function is obtained by remov- 
ing one row and column from the determinantal matrix. 
G a (iuj n ) simply relates to M = A -1 by 12 - 15 

G imp (iujn) =-\Y. M^e^ niTi ~ Ts) (4) 

hj 

Alternatively, one can simulate the impurity Green's 
function from the cluster trace at each Monte Carlo 
step; 15 that is, 



Gimp (j'^n) 



1 



^{ C {T)c^)dT 



fk 11 dTelU1 " T EW'I^M^c^) (5) 



Here, \cf>) is the eigenstate of the Anderson impurity 
model, in terms of which the full partition function can 
be written as Z — J2<p e ~ ■ For each specific config- 
uration Ck sampled in the CT-HYB, this expression has 
the following form: 



x 



Y,{m\e-^T[c{r)T^\ m ) (6) 



The explicit form of the determinant is given in Eqn. (3). 
T[ and Th 1 are the left and right lists of cluster operators 
c(t), respectively, with the constraint t; + i < t < t\. The 
partition function corresponding to the configuration Ck 
is given as 

Z Ck = Z loc Z b ^2{m\e-^Tl x Tt +1 \m)Det Ck 

m 

= Z loc Z h Tr{C k )Det c K (7) 
By combining the above two equations, we have 
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with the notation = (jn\T{\n) and 

e (iui n +E n -E p )Ti e (iu n + E n -E p )n + 1 



The ratio Zc k /Z is the probability of configuration Ck 
being sampled in the Monte Carlo simulation. 

When k a is small, we measure Gi mp directly from the 
cluster trace, 15 , that is, Eq. (8). Although this scheme 
is not very fast, it is more stable than Eq. (4). When k a 
is large and Eq. (4) is used in the simulation, the high- 
frequency parts of Gi m p converge much slower and con- 
tains more statistical errors than the low-frequency parts. 
As a result, the corresponding self-energy can be fluctu- 
ating at large w„ . As already pointed out in the Introduc- 
tion, the correct high-frequency behavior of Y>i mp (iuj n ) is 
crucial for the CT-HYB. Thus, special attention has to 
be paid to get rid of the noises in the self-energy data. 

To the best of our knowledge, three schemes have been 
proposed for dealing with this problem. (1) Noise fil- 
tering. One can either smooth the noises at r ~ j3/2 
by averaging Gi mp (r) over a small range of r (see Refs. 
11 and 12) or apply the orthogonal polynomial filter- 
ing routine recently proposed by Boehnke et al. le to 
achieve a smooth C?i mp (r) for all r e [0, /?). By care- 
fully choosing the order of the orthogonal polynomials, 
the impurity self-energy becomes smooth for all Matsub- 
ara frequencies. (2) Replacing the high-frequency tail 
of Sj mp (io; n ) with some well-behaving function. This 
function can be either the self-energy, calculated from 
a weak-coupling perturbation expansion, or the moment 
expansion of the Green's function. 17,18 Such a replace- 
ment provides a smoothly behaving high-frequency tail 
of the self-energy function. However, the correspond- 
ing expression usually becomes complicated in the multi- 
orbital case and relies on the explicit form of the inter- 
action term. (3) Measuring Gi m p(r) from higher order 
correlation functions. 19 This method becomes advanta- 
geous for the density-density type interaction, for which 
the "segment picture" 11 can be used. For general type 
interactions, numerical cost has to be paid to calculate 
additional correlators. 

Here, we propose a simple and stable scheme which 
does not rely on any direct noise filtering of Gi mp {r) and 
does not introduce any numerical cost to the CT-HYB 
simulations. This method docs not depend on the explicit 
form of the interaction term and remains efficient in the 
multiorbital calculations. The basic idea is to determine 
an approximate self-energy function by performing the 
perturbation expansion around the atomic limit, using 
the 'dual-transformation'. As we will see later on, such a 
method generates systematic improvements to the atomic 
self-energy. The first-order expansion term already gives 
considerable corrections and reproduces the correct high- 
frequency behavior of T, imp (ito n ) . 

The expansion around the atomic limit has been stud- 
ied before. 20 In the strong-coupling region, this method 
yields results comparable to the numerical exact QMC 



results. Here, we use an elegant and different way, that 
is, the "dual transformation" . 21 This transformation has 
been used in the construction of the dual-fermion (DF) 
method, which gives an action well behaving in both 
the weak- and the strong-coupling limits. Thus, our 
perturbation expansion actually also works in the weak- 
coupling region. 

The impurity model has the following action: 

S[c*,c] = Sjmp[c*, c] + 7^ y, c* A a (iu) n )c a (10) 

n a 

In the "dual transformation" , new variables /* , / are in- 
troduced to rewrite the hybridization term in the follow- 
ing way: 

e <A a ( iw „) Co det^p 1 
or 

= j e -«(<f*+fZ c «)-±£^)fiu V [f*j}. (11) 

The complex number a can be arbitrary in the above ex- 
pression. In Ref. 21, it is taken as the impurity Green's 
function. This makes the correlator of the dual variables, 
that is, G d = —(fa fa), behaves like the one-particle 
Green's function, which decreases as l/iiv n for large u n . 
For simplicity, we take a as one. Although in this case, 
the dual variables can not be interpreted as fermions, the 
impurity Green's function remains the same. 

Integrating out the c variable, the full action becomes 
a functional which only depends on variables /*, /, that 
is, 

Z = Z f Z b J V[r,f}e-^K A ^f° 

= Z f Z b J X>[/*, /]e~ S» /; G 5'~ l /--vi* } A^/a/X( 12 ) 

where Gq is given as [G"' — A" 1 ] -1 . The effective inter- 
action of dual variables turns out to be the reducible 
four-point correlations of the atomic system, that is, 
Vt ] = X& i si-Si,20G&G& + 6 1A J3GftG%, with G at be- 
ing the atomic Green's function. 

Since the dual transformation is mathematically ex- 
act, the two different actions which depend on only c 
variables [i.e., Eq. (10)] and /variables [i.e., Eq. (12)] are 
equivalent. Thus, we can obtain an exact relation be- 
tween the correlators G a and Gf from differentiating the 
two actions with respect to A a . This yields: 

G„ = -A- 1 -A- 1 G^A- 1 , (13) 

where G„ is obtained from the Dyson equation, that is, 
G d a = [G{f - Sf]" 1 . is the self-energy function of the 
dual variables. The expression of Xi2^4 can be found in 
the literature, (e.g., Refs. 22-24). If the interaction of the 
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dual variables in Eq. (12) is neglected, the atomic self- 
energy will be recovered. This can be seen by inserting 
into Eq. (13). We have 

G a (iu n )=G?/(l-A a G: t ). (14) 

Then, from the Dyson equation, we immediately see that 

J^* = iu n + li-A a -G- 1 (iu n ) 

= uj n + li-l/G? = £* (15) 

Thus, one can imagine the interaction term in Eq. (12) 
will generate systematic corrections to the atomic self- 
energy. 

By including the interaction and further restricting the 
calculation of to the first order, we have 

K^n) = -i E E K«S(^; iu' n )G d b «) (16) 
P b ui' n 

In this equation, only the element V^'.^ £12^34 is re ~ 
quired. Additionally, this calculation can be further ac- 
celerated by employing the look-up routine and the sym- 
metry of Xi2;34, which is shown in Sec. II B. By doing 
so, the perturbation expansion remains very efficienct in 
multi-orbital calculations. 

As a benchmark, we first apply the dual expansion 
scheme by restudying the Bethe lattice with different 
bandwidths, that is, W<z = 2Wi, where the orbital- 
selective Mott transition can happen. 11,25-35 We directly 
solved the DMFT equation with the high-frequency sup- 
plemented self-energy function, instead of using Eq. (20) 
in Ref. 12. Our self-energy data in Fig. 1 is identical to 
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FIG. 1. Benchmark: imaginary part of the impurity self- 
energies for fit\ — 50, U/ti = 4, J/U — 0.25 in unit of t\. 
Both the narrow and the wide bands are metallic. The dual 
expansion gives two different asymptotic behaviors of the self- 
energy for two different bands, as expected. However, the 
atomic self-energy does not have such a resolution. 

those in Fig. 12 of Ref. 12, meaning that the dual ex- 
pansion method is reliable to produce the high-frequency 
tail of the self-energy and can be used in the CT-HYB 



for solving impurity problems. To see the performance of 
the dual expansion method for a finite spatial-dimension 
problem, in Fig. 2 we show the comparison of the self- 
energy function calculated for a two-orbital Hubbard 
model in two dimension [see the Hamiltonian in Eq. (1)]. 
The improvement from the dual expansion is clearly seen 
from the agreement between the CT-HYB and the dual 
expansion results. Increasing the hybridization strength, 
this agreement becomes even better. Thus, a smaller 
number of Matsubara frequencies is required to simu- 
late in such a case. However, the atomic self-energy has 
a larger deviation from the CT-HYB results for smaller 
u) n . Similar ideas were used to formulate effective impu- 
rity solvers 22 ' 24 for the DMFT. We use it here to get the 
correct high-frequency tail of the impurity self-energy, 
while still keeping the low-frequency self-energy function 
simulated from the QMC. This method only needs the 
hybridization function at each DMFT iteration. The 
dual-expansion can be carried out independently of the 
CT-HYB simulation. Thus, it does not introduce addi- 
tional numerical cost to the CT-HYB, which is another 
essential difference with respect to previous works. 15-19 




2 4 6 8 10 

FIG. 2. The comparison of the impurity self-energy calculated 
from the CT-HYB, atomic Hubbard model and the dual ex- 
pansion method. The parameter sets for the two-orbital Hub- 
bard model are /3t = 50,U/t = 4.0, J/U = 0.25, V/t = 1.0. 
See the text for more details. 



B. Four-point correlation function Xi2 ; 34 

The dual expansion, discussed in the above section, 
requires the knowledge of the atomic four-point corre- 
lation function XiI-34- I n ^ ne multiorbital case, such a 
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calculation can be hard since the large-dimensional ma- 
trix multiplication is time consuming. In this case, one 
can again use the block diagonal form of the Hamiltonian 
matrix and employ the look-up routine as we did in the 
trace calculation. Here, we want to further simplify the 
calculation by employing the symmetry of Xi2-34- 

Such 

a symmetry turns out to be also very useful in the sim- 
ulation of the impurity four-point correlation function 
Xi2-34- Thus, in this section we try to keep our discus- 
sion general. We start from the simulation of the x™34 
in the CT-HYB and discuss the symmetries of it after- 
ward. The same symmetry requirements are satisfied by 

X?2;34 as welL 

Although in the CT-HYB, Wick's theorem apparently 

is not supported by the impurity action, the four-point 

correlation function can be simulated by removing two 

rows and two columns from the determinant matrix, 

which results in an expression analogous to those for 

the CT-INT and the CT-AUX. Effectively, one can still 

simulate the four-point correlation function as if Wick's 

theorem holds. Here, we use the following notation to 

symbolically represent this expression: 

Xl2;34 : ~ \ C 1 C 2 C 3 C 4/ 
= ( C 1 C 2>( C 3 C 4> - ( C 1 C 4)( C 3 C 2) 

= ffi2 H > w 2 )S 3 4 ( w 3 » u A ) - g u H , w 4 )g 32 (cj 3 , W 2 ) (17) 

where labels 12; 34 represent "orbitals, sites, spins," etc. 
In the CT-HYB, the two-frequency dependent propaga- 
tors g a p(u>,w') is given as 

Saj9 (a; 1>W2 ) = -I^ e -^M<f e"^. (18) 

id 

It has the following symmetry in Matsubara frequency 
space: 

g a p(ui,U2) = g* aP (-uJi,-uj2), (19) 

which reduces the numerical effort by a factor of two. A 
similar symmetry is also satisfied by X- 

Xl2;34 = X?2;34( W ' W = Xl"it{-V, -U>'). (20) 

In what follows, we denote u> = uji, u> + fl = u>2, + fl = 
0J3, a/ — UI4. Equation (20) says, only for fl > 0, x needs 
to be simulated. 

Symmetry (20) relates the positive frequencies to the 
corresponding negative frequencies of X- It is also possi- 
ble to find symmetries which connect different u>, a/ in 
the same fl sector. This can be achieved via the fact 
that Xi2;34 is antisymmetric under the exchange 1 3 
and 2 <^> 4: 

X34;12(^' + fl,uj + fl; -ft) = Xi2:34(w,u/;f2). (21) 
Combining Eq. (21) with Eq. (20), we have 

X34;12(-^' - ft, -U - fl; fl) = X l 2 - M {"> ^ (22) 



Given the spin configurations of different x channels, we 
find Xi2-34 <T satisfies both symmetries in Eqs. (20) and 
(22). However, Xi2-34 om y satisfies the symmetry shown 
in Eq. (20) and the following relation: 

Xl2;34 — X34;12 • \ Z6 ) 

One can implement the symmetries in Eqs. (20) and 
(22) as follows. (1) x i t^^-W, and x are simu- 
lated only for fl > 0. (2) For each specific fl considered, 
Eq. (22) is further applied to x^'^ and x^'^- Only for 
parts of the frequency points in this fi-sector do x^'^ 
and x^'^ need to be simulated. (3) At the end of the 
calculation, fl < components are calculated through 
Eq. (21). (4) x U '' n is calculated by Eq. (23). In ad- 
dition to the symmetries shown in Eqs. (20) and (22), 
it is possible to find more symmetries to relate different 
frequency sectors. 

Before finishing this section, we want to note that 
the four-point correlation function is useful not only for 
the physical response function and the dual-expansion 
scheme, but also relates closely with the extension of the 
DMFT. In the DF method 21 and the dynamical vertex 
approximation (DrA), 23 the nonlocal self-energy is con- 
structed from the impurity two-particle vertices. 



III. APPLICATION 

As a typical application, we consider here a two-orbital 
Hubbard model [see the Hamiltonian in Eq. (1)], with ro- 
tationally invariant interactions, i.e. U' = U — 2 J, U" = 
U' — J. To make a link with realistic material systems, 
this multiorbital Hubbard model can be viewed as an ef- 
fective model for the e g -orbital systems. The rotational 
invariance of the interaction term is not obligatory in 
the CT-HYB solver; here, we use it only as one possible 
situation. By making use of the DMFT, the two-orbital 
Hubbard model has been studied by many groups. 12,3 
These calculations are either based on a semicircular den- 
sity of states, which corresponds to the Bethe lattice, or 
they employ an impurity solver with certain limitations 
in temperature or interaction strength. Here, we solve 
the DMFT equation at finite dimension and tempera- 
tures. In these cases, the DMFT loop cannot be closed 
by a simple relation in the imaginary-time space like on 
the Bethe lattice. Thus, our dual-expansion method dis- 
cussed in Sec. II A turns to be a decisive tool. Our calcu- 
lations are mainly performed on ordinary desktop com- 
puters. 

Compared to the single-orbital case, two issues in a 
multiorbital model are of obvious interests: 

(1) What is the effect of the orbital fluctuations ? The 
general believe is, that it is competitive to the Coulomb 
interaction. As a result, the metallic state can be stabi- 
lized up to a large interaction value 31,44 . 

(2) How does the Hund's coupling modify the transition 
from the metal to Mott insulator (MIT)? It is known 
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FIG. 3. The phase diagram of the two-orbital Hubbard model 
at half filling. The MIT at fit = 50 for two values of J/U 
are shown as histograms. The two local density of states 
on the right-hand side correspond to the two solutions for 
U/t = 8.4, J/U = 0.1. 
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FIG. 4. Behavior of the impurity self-energy and Green's 
function around the metal-to-band-insulator transition as 
functions of the hybridization strength V/t. 



that the two-orbital Hubbard model behaves quite dif- 
ferently with and without J. 36 ' 37 

The phase diagrams of the two-orbital Hubbard model 
can be found in Refs. 37 and 42. Here we study, in par- 
ticular, the coexistence region for different values of J in 
Fig. 3 (a), which indicates the MIT is of first order. Com- 
pared to the phase diagrams for the Bethe lattice, 37 ' 42 the 
reduction of the spatial dimension does not change signif- 
icantly the critical Coulomb interaction value of U c when 
it is normalized by the full bandwidth. However, U c be- 
comes larger compared to the single-orbital model, which 
confirms that the orbital fluctuation stabilize the metal- 
lic phase. With the increase of the Hund's rule coupling 
J, we found the coexistence region to become smaller. 
For the two values of J/U in our calculations, the reduc- 
tion is about 0.2 eV. On the other hand, Bulla et al. 45 
found, for J/U > 0.25, the transition to be of second 
order. At J/U = 0.25, our results show that the coexis- 
tence region still has a reasonably large width. Thus, we 
believe that even for J/U > 0.25, the MIT remains first 
order. Whether, the coexistence region completely dis- 
appears with the further increase of J/U deserves more 
investigations. 

On the right-hand side of Fig. 3, two different solutions 
of the local density of state, that is, A(oj), are displayed 
for U/t — 8.4. They correspond to the metallic, see Fig. 
3 (b), and insulating states [see Fig. 3 (c)] in the coex- 
istence region. A(u) is obtained by using the stochastic 
analytical continuation directly on the Matsubara data 
of G imp (iu! n ). i6 

In Fig. 4, the typical behavior of the metal-to-band- 
insulator transition is shown by calculating the impu- 
rity Green's function and the corresponding self-energy 
as a function of the hybridization. Increasing the hy- 
bridization V/t tends to open a band gap. Furthermore, 
with the increase of V/t, the impurity Green's function 



at the lowest Matsubara frequency becomes smaller and 
finally approaches zero [see Fig. 4 (a)]. The metal- 
to-band-insulator transition happens somewhere between 
V/t = 2.5 and 3.0 for U/t = 4. This transition is not vis- 
ible from the self-energy plot, where Y>i mp (iu! n ) behaves 
similarly for different values of V/t. The slope, that is, 
dY*i mp (u})/duj\ U j , remains negative for all hybridization 
strengths [see Fig. 4 (b)]. In contrast, the slope of the 
local Green's function around loq has different signs be- 
fore and after the metal-insulator (band) transition. 

Increasing further the value of U/t strengthens both 
the intra- and intcrorbital interactions. Finally, for val- 



-0.02 
3" -0.04 
£ -0.06 

-0.08 



— • V/t =0.0 

t — • V/t =0.125 

»— » V/t =0.25 

•— V/t =0.4 

■— ■ V/t =0.6 

° ° V/t =0.75 





FIG. 5. Similar to Fig. 4, but with different interaction 
strength U/t = 9, where a Mott-insulating state is found 
at A/i = 0,0.125,0.25. The increase of the hybridization 
between two orbitals greatly changes the behavior of the self- 
energy, while it leaves the one-particle Green's function es- 
sentially unchanged. 
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ues of U/t of the order of the noninteracting bandwidth, 
the metal to band-insulator transition is replaced by the 
Mott-insulator-to-band-insulator crossover as a function 
of the hybridization strength A/t. This behavior is dis- 
played in Fig. 5. In contrast to the metal to band- 
insulator transition shown in Fig. 4, in Fig. 5 (with a 
choice of U/t = 9), the local Green's function stays nearly 
unchanged under modifying the hybridization strength 
V/t, that is, Gi m p shows an insulating behavior for all 
values of V/t. However, for different values of V/t, the 
insulating nature is indeed different. This can be seen 
from the variation of the self-energy function shown in 
the right-hand side of Fig. 5. Increasing V/t, results 
in increasing of <9£(w)/9w| Wo for any finite V/t, indicat- 
ing the crossover from Mott-insulator to band-insulator 
behavior. 47 

By applying the symmetries presented in Sec. II B, 
we show the results for the interorbital and intraorbital 
reducible spin susceptibilities in Fig. 6 for fit = 20, U/t = 
6, J/U = 0.25, and V/t — 0, with a, b the orbital indices: 

Xn (w»,wj = -{X ab ,n ~ X a b,n ) ( 24 ) 

X^ m ' ab (uj n ,uj' n ) are the impurity susceptibilities with 
the subtraction of the impurity bubble susceptibilities. 
They are plotted as functions of the two fermionic fre- 
quencies co n ,uj' n for fixed O = 0. While here only the 
51 = component is given, the implementation discussed 
in Sec. II B works for any value of f2. Figures 6(a) 
and 6(b) refer to the three-dimensional (3D) plots of 
XQ ln ' ab (oJ n , ^4)' the corresponding 2D top- view plots are 
shown in Figs. 6(c) and 6(d). Based on the CT-HYB, the 
four-point correlation functions were recently also calcu- 
lated for the effective one and four-orbital systems 16,48 
for different problems. Another efficient and stable, but 
approximate, algorithm can be found in Ref. 49. 

From Fig. 6, we see that the reducible two-particle sus- 
ceptibility XQ^' ab (uJ n , u> n ') decays rather fast as a func- 
tion of u> n and uj n i . The dominant contribution comes 
from the elements with u) n = 0, or uj n > = 0, or w„ = u) n i. 
For our parameter set, the interorbital spin susceptibil- 
ity shows a sharper structure than the intraorbital one, 
which can be viewed as a precursor of the possible orbital 
antiferromagnetic order. 



IV. CONCLUSION 

In this paper, we showed how the high-frequency tail of 
the self-energy can be calculated in a controlled manner 
from the dual transformation in CT-HYB. This scheme 
provides an efficient recipe for finite-dimension DMFT 
studies when taking the CT-HYB as an impurity solver. 
Our procedure is based on a Matsubara frequency space 
simulation and produces more moments from the dual ex- 
pansion. Thus, it generates an improved high-frequency 
self-energy behavior. Most importantly, it does not intro- 




FIG. 6. The interorbital, that is, Xf2=o> an d intraorbital, that 
is, Xn=oi components of the reducible impurity two-particle 
susceptibility for fit = 20, U/t = 6, J/U = 0.25. 



duce any additional numerical cost to the runtime simu- 
lation. We also simulated the four-point correlation func- 
tion for different spin configurations in the particle-hole 
channel. To this end, we implemented different symme- 
tries to reduce the memory and CPU requirements with- 
out losing accuracy. 

As a first application, we demonstrated the usefulness 
of our method for a two-orbital model with a general on- 
site interaction. From this study, we deduced a substan- 
tial influence of the Hund's rule coupling on the metal- 
insulator transition phase diagram, especially on the co- 
existence region. In particular, we find that for any finite 
value of J/t, the MIT stays first order. 

Our scheme is also of particular use for connecting the 
DF method, which many be viewed as a nonlocal exten- 
sion of the DMFT, with a priori DFT techniques. A 
multiorbital DF calculation will be especially interesting 
and rewarding for the DFT + DF study of material sys- 
tems. In such study, the CT-HYB effectively works on 
an impurity problem with the DFT dispersions as input. 
Thus, one has a good control on the "minus-sign" prob- 
lem. The high-momentum resolution, provided by the 
DF algorithm, makes the result ready to be compared 
with experiments, such as ARPES data. 
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